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ABSTRACT 

We report the first results of an ^yVM-Newton monitoring campaign of the NGC 6231 
open cluster in the Sco OB 1 association. This first paper focuses on the massive col- 
liding wind binary HD 152248, which is the brightest X-ray source of the cluster. The 
campaign, with a total duration of ISOksec, was split into six separate observations, 
following the orbital motion of HD 152248. The X-ray flux from this system presents 
a clear, asymmetric modulation with the phase and ranges from 0.73 to 1.18 10~^^ erg 
s~^ cm~^ in the [0.5-10.0 keV] energy band. The maximum of the emission is reached 
slightly after apastron. The EPIC spectra are quite soft and peak around 0.8-0.9 keV. 
We characterize their shape using several combinations of MEKAL models and power- 
law spectra and we detect significant spectral variability in the [0.5-2.5 keV] energy 
band. 

We also perform 2-D hydrodynamical simulations using different sets of param- 
eters that closely reproduce the physical and orbital configuration of the HD 152248 
system at the time of the six ^VWl-Newton pointings. This allows a direct confronta- 
tion of the model predictions with the constraints deduced from the X-ray observations 
of the system. We show that the observed variation of the flux can be explained by 
a variation of the X-ray emission from the colliding wind zone, diluted by the softer 
X-ray contribution of the two 0-type stars of the system. Our simulations also reveal 
that the interaction region of HD 152248 should be highly unstable, giving rise to shells 
of dense gas that are separated by low density regions. 

Finally, we perform a search for short-term variability in the light curves of the 
system and we show that trends are present within several of the 30 ksec exposures of 
our campaign. Further, most of these trends are in good agreement with the orbital 
motion and provide a direct constraint on the first order derivative of the flux. In 
the same context, we also search for long-range correlations in the X-ray data of the 
system, but we only marginally detect them in the high energy tail of the signal. 

Key words: stars: individual: HD 152248 - binaries: close - stars: early-type - stars: 
winds - X-rays: individual: HD 152248 - X-rays: stars. 



1 INTRODUCTION 

Early-type stars of spectral type O and their evolved 
descendants, the Wolf-Rayet stars, are among the most 
luminous and hottest objects of the Milky Way. Charac- 
terized by their strong and powerful winds, and their huge 
mass loss rates, they are also known to be X-ray emitters 
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since the late 70's and the advent of the X-ray observatories 
equipped with Wolter focusing mirrors. Though with some 
dispersion, the X-ray emission from single O and early B 
stars approximately scales with their bolometric luminosity: 
Lx ~ 10~''Lboi (Berghofer et al. 1997) and displays a 
soft thermal spectrum {kT ~ 0.2-1.0 keV). The exact 
mechanism of this X-ray emission, however, is still not fully 
understood but it is generally thought that, in the lower 
layers of the winds, the plasma is heated in shocks that 
arise from small-scale wind structures that grow out from 
line-driven wind instabilities (Feldmeier, Puis & Pauldrach 
1997; Dessart & Owocki 2003). 
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Beyond this general scheme, it is now established that 
massive binaries often display an additional X-ray lumi- 
nosity compared to single stars of the same spectral types 
and luminosity cleisses (Chlebowski & Garmany 1991). The 
extra X-ray component is usually attributed to emission 
from a hot plasma that results from the hydrodynamical 
collision of the winds of the two stars. The standard collision 
model predicts the interaction region to be limited by two 
curved shock-surfaces between which the gas is heated up 
to temperatures of a few W'^K. From the observational 
point of view, this phenomenon can manifest itself i) as a 
modulation of the X-ray emission due to the variation of the 
line-of-sight opacity with the orbital motion, ii) in an eccen- 
tric binary, as a variation of the intrinsic emission due to a 
variation of the separation between the two components, as 
the densities of the colliding material and, in close binaries, 
the pre-shock wind velocities depend on the distance to the 
originating star. Though the X-ray domain is probably the 
most adequate to study the wind interaction region in close 
binary systems, the observational signature of the collision 
is not restricted to high energies. Prom the UV and optical 
domains for example, the interaction may be traced in 
emission lines that arc produced by recombination of the 
cooling gas that escapes the interaction region after the 
shocks (see e.g. Sana, Rauw & Gosset 2001; Rauw et al. 
2002b). 

In addition to the improvements in observational 
capabilities, another tool for probing the physics of wind- 
wind collisions has undergone important development in 
the past 10 to 15 years. Namely Computational Fluid 
Dynamic (CFD) applied to the hydrodynamical problem 
of wind-wind collision holds out promise of insight into the 
physics of the phenomenon (Stevens, Blondin & Pollock 
1992; Pittard & Stevens 1997; Henley, Stevens & Pittard 
2003; De Becker et al. 2004a). For example, from the work 
of Stevens et al. (1992, SBP92 hereafter), it is obvious that 
the shocks and the interaction zone properties arc strongly 
dependent on the characteristics of the fluids upstream 
the shock. In a first step to classify the shocks, SBP92 
introduced a parameter x = icooi/ifiow that represents 
the ratio of the characteristic cooling time of the fluid 
downstream the shock (icooi) to the characteristic flow 
time (tflow)- Based on this simple criterion, two extreme 
behaviours might be distinguished. On one hand the gas 
escapes the interaction region with almost no significant 
cooling (x >> 1) and the phenomenon is mostly adiabatic. 
On the other hand, whenever x ^ li the radiative cooling 
then plays a crucial role. In the adiabatic case, the escaping 
flow is expected to be relatively smooth and the region 
in botwccu the shocks is quite large. When the flow is 
dominated by radiative cooling, however, the shock region 
collapses and the steady state is disrupted by Kelvin- 
Helmholtz and thermal instabilities (SBP92). 

However until recently very few X-ray observations of 
O-l-O Colliding Wind Binaries (CWBs) had a sufficient 
quality (counts, time resolution, ...) to allow a detailed 
comparison with predictions from CFD calculations. With 
the advent of the XMM-Newton and Chandra X-ray 
observatories, the hot star community was offered for the 
first time the possibility of acquiring high quality data 



that would provide a strong tost of the latest models. This 
should therefore noticeably contribute to improving our 
understanding of X-ray emission from early-type stars, and 
more specifically of the winds of these particular objects. 
However, despite their interest, still very few 0+0 binaries 
have been the target of such an X-ray monitoring to detect 
the expected phase-locked behaviour of the emission. 

HD 93403 is an 05.51 -I- 07V binary with an orbital pe- 
riod close to 15.1 days and an eccentricity of 0.234 (Rauw 
et al. 2000). An XMM-Newton campaign revealed a phcise 
locked behaviour of the X-ray flux that is mostly consistent 
with a l/Dscp dependence, where Dsop is the separation be- 
tween the two stars, as expected for an adiabatic wind-wind 
interaction (Rauw et al. 2002b). 

A similar behaviour is observed in the case of HD 93205, 
an eccentric system (e = 0.37, Porb = 6.08 days) consisting 
of an 03 V primary and an 08 V secondary, in the Tr 16 clus- 
ter near rj Car. Antokhin et al. (2004) analyze five XMM- 
Newton observations of the region around rj Car. They re- 
port a clear phase-locked modulation of the flux, with maxi- 
mum emission occurring at periastron and minimum around 
apastron. 

In the case of the 07V-I-07V binary HD 159176 
(Porb = 3.367 days), De Becker et al. (2004a) obtained a 
single observation with XMM-Newton. The system was 
found to display a rather modest X-ray overluminosity 
attributed to a wind collision. The authors further find that 
only part of the wind kinetic power is actually emitted in 
the X-ray domain as a result of the wind interaction. 

In this context, we undertook a phase-resolved 
XMM-Newton campaign targeting the early-type binary 
HD 152248. This system lies at the centre of the young open 
cluster NGC 6231 in the core of the Sco OB 1 association. 
It is a close SB2 eclipsing binary consisting of two almost 
identical O-type giants with a period close to 5.8 days. 
The orbital and physical properties of the system were 
established by Sana et al. (2001, Paper I) from medium and 
high resolution optical spectroscopy. In Paper I, we showed 
that HD 152248 is an 07.5(f) III + 07(f) III system with 
the 07.5 HI primary being slightly less massive than the 
secondary. The convention that refers to the less massive 
component as the primary star historically results from the 
light curve of the system that displays a deeper minimum 
during the occultation of the 07.5 star. Table 1 provides 
the physical and orbital parameters of the system. Though 
HD 152248 could not be resolved with the EINSTEIN 
observatory due to the limited spatial resolution of the 
satellite and the crowded nature of the field, HD 152248 
was observed by the Rontgen Satellite (ROSAT) and clear 
variations of the flux with the orbital phase could be iden- 
tified (Corcoran 1996). Paper I also presents strong clues 
that HD 152248 harbours a colliding wind phenomenon 
that could be tracked in the optical domain through the 
phase-locked profile variations of the He II A 4686 and Ha 
emissions. 

It therefore appeared extremely promising to compare 
the latest X-ray observations of this object with the 
predictions of recent hydrodynamical models. Indeed the 
quality of the HD 152248 XMM-Newton data combined 
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Table 1. Orbital and physical parameters of the HD 152248 bi- 
nary. The usual notations have been adopted. To is the time of 
periastron passage and is adopted as phase cp = 0.0. Most of the 
quoted values are from Paper I. Whenever this is not the case, 
reference to the original paper is provided. 
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with the very good knowledge of the system and the high 
level of constraints on its physical and geometrical param- 
eters (including the usually rather tricky determination 
of the incUnation) provide the opportunity i) to perform 
numerical simulations corresponding to a well constrained 
configuration, and consequently very representative of 
the HD 152248 geometry, ii) to provide predictions that 
are suitable for a direct confrontation with observational 
results, iii) to test the state-of-the-art CWB models, iv) 
to gain insight into the physics of wind-wind collisions 
and v) to uncover new routes to improve hydrodynamical 
simulations of close binary systems. In this paper, we thus 
discuss the X-ray emission from the HD 152248 system 
both from the observational point of view and by means 
of hydrodynamical calculations of wind-wind collisions. 
The analysis of the other X-ray sources of the field will 
be addressed in subsequent papers. We further focus on 
the EPIC observations of the system and on a detailed 
comparison with corresponding numerical simulations. We 
do not tackle here the analysis of the RGS data. Indeed, as 
the field around HD 152248 is relatively crowded, the RGS 
spectra are contaminated by rather bright neighbouring 
sources and a simultaneous analysis of the EPIC data of 
the contaminants as well as specific techniques to estimate 
the level of the contamination are required. This is beyond 
the scope of this study and will be deferred to a future work. 

The present paper is organised as follows. The next sec- 
tion presents a description of the XMM-Newton campaign 



Figure 1. Scale sketch, in the orbital plane of the system, rep- 
resenting the configuration of the HD 152248 binary at the time 
of the six XMM-Newton pointings. The primary star is in dark 
grey while the secondary is represented in light grey. Arrows at 
Icfthand indicate the projection, onto the orbital plane, of the line 
of sight of the satellite towards the system. 



Figure 2. MOSl (upper panel) and pn (lower panel) images of 
the central part of the NGC6231 cluster in the energy band [0.5- 
2.0 keV]. Adopted extraction regions for HD 152248 and for the 
background are also indicated. 



and the data reduction. Sect. 3 is focused on XMM EPIC 
spectra analysis and, in Sect. 4, we describe the performed 
hydrodynamical simulations. We also provide a detailed 
comparison between observation and simulation results. 
The next section (Sect. 5) investigates the short-term 
variability of the observed X-ray emission. Finally, Sect. 6 
presents a summary of the main results of this paper as 
well as the conclusions of the present work. 

2 OBSERVATIONS AND DATA REDUCTION 

2.1 The XM.M-Newton observing campaign 

The X-ray observational campaign was designed in the 
framework of the X.MM-Newton guaranteed time of obser- 
vation of the XMM-OM consortium. It initially consisted 
of six SOksec pointings towards the young open cluster 
NGC6231 that shelters the HD 152248 binary. These six 
observations were spread over the 5.816 day period of the 
system (see Fig. 1) to monitor the expected phase-locked 
behaviour of the emission. In September 2001, the XMM- 
Newton satellite successfully performed the six observations 
within satellite revolutions 319 to 321, i.e. within a single 
orbital cycle of HD 152248. Two of the six data sets were 
unfortunately affected by soft proton fiares which effectively 
reduced the observing time by about one third. Table 2 
gives the journal of the X-ray campaign dedicated to the 
cluster NGC6231. All three EPIC instruments were oper- 
ated in the Full Frame mode together with the thick filter 
to reject optical light. The RGS spectrographs were run in 
the Standard Spectroscopic mode. Due to the brightness 
of the objects in the field of view (FOV) of NGC6231, the 
Optical Monitor was switched off throughout the campaign. 



2.2 XMM EPIC data reduction 

The EPIC Observation Data Files (ODFs) were processed 
using the XMM-Science Analysis System (SAS) v5.2 
implemented on our computers in Liege. We applied the 
emproc and epproc pipeline chains respectively to the MOS 
and pn raw data to generate proper event list files. No 
indication of pile-up was found in the data. We then only 
considered events with patterns 0-12 (resp. 0) for MOS 
(resp. pn) instruments and we applied the filtering criteria 
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Table 2. Journal of the XMM-Newton observations of HD 152248. Columns 2 and 3 give the spacecraft revolution number and the 
observation ID. The mean Julian Day (JD) is reported in Col. 4. The next three columns list the performed exposure times for the EPIC 
MOS, EPIC pn and RGS instruments. The last column provides the orbital phase of HD 152248 for each XMM-Newton observation at 
mid-exposure, according to the ephemerides given in Table 1. The quoted uncertainties stand for the phase intervals that correspond to 
the durations of each observation. 
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28.5 


31.6 
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6 


321 


0109490601 


2162.726 


32.9 


30.3 


33.5 


0.31 ± 0.03 



XMMEA_EM and XMMEA_EP as recommended by the 
Science Operation Centre (SOC) technical note XMM-PS- 
TN-43v2.0. For each pointing, we rejected periods affected 
by soft proton flares. For this purpose, we built light curves 
at energies above lOkoV^ and discarded high background 
observing periods on the basis of an empirically derived 
threshold. The so-dcfiucd GTIs (Good Time Intervals) were 
used to produce adequate X-ray event lists from which we 
extracted images and spectra. 

For the purpose of spectral analysis, we adopted a circu- 
lar extraction region centered on HD 152248. However, the 
NGC 6231 cluster around this object is relatively crowded 
(see Fig. 2) and we could not adopt an extraction radius 
larger than about 40", which corresponds to a fractional en- 
circled energy of about 85 per cent. This extraction region 
was used for the three EPIC instruments throughout the 
data treatment of the observing campaign. For the same rear 
son, we could not adopt an annular region around the source 
to estimate the background level. After several tests, we se- 
lected a circular region located outside the crowded part of 
the cluster (see Fig. 2). As a complementary test, we also de- 
termined the background level using the blanksky files cre- 
ated from deep field observations. These arc described in the 
SOC technical note CAL-TN-0016-1-1 v2.0. However these 
observations were performed with the thin filter whereas 
our data were obtained with the thick one. As the back- 
ground sky level may be dependent on this factor, we tried 
to account for it: we extracted the thin background spec- 
tra from the blanksky files using a region identical to the 
extraction area adopted for the HD 152248 binary. We then 
adjusted an empirical model to the thin background and 
generated an equivalent TinCK background using the fake 
XSPEC command. For this operation, we used the arf+rmf 
response file provided by the SOC. This latter approach was 
limited to the EPIC MOS data sets as a proper correspond- 
ing pn BLANKSKY file was unavailable. However, the two 
background techniques give consistent results, well within 
the mutual error bars. This supports the a priori idea that, 
due to the brightness of HD 152248, the background deter- 
mination is not a critical issue for this object. In the fol- 



^ Expressed in Pulse Invariant (PI) channel numbers and consid- 
ering that 1 PI channel approximately corresponds to 1 eV, the 
adopted criterion is actually PI>10 000. 



Figure 3. EPIC MOS and pn light curves (in cnt s'^) of the 
HD 152248 system in the 0.5-10.0 keV energy band. Horizontal 
bars represent the effective duration of the exposure while vertical 
ones axe the 1-cr uncertainties on the count rates. 

lowing, we adopt the background as determined from our 
observations. We generated adequate rm/ and ar/ files using 
the rmfgen and arfgen commands. We also used the ma- 
trices provided by the SOC. Again no difference was found 
between the spectra obtained in both ways. Finally, back- 
ground corrected spectra were produced using the grppha 
command of the ftools package. 



3 XMM EPIC OBSERVATIONS 
3.1 Light curves 

As a major aim of this project is to monitor the X-ray 
variability of the HD 152248 system, we first extracted 
broad band light curves from the EPIC data. Plots of the 
count rates against the phase (see Fig. 3) clearly reveal an 
important enhancement of the global emission by about 60 
per cent at the apastron passage. The maximum is reached 
at = 0.66, slightly after apastron {<f) = 0.51), and is fol- 
lowed by a rapid decrease that slows down after periastron 
passage. To refine our analysis, we selected four energy 
bands: very soft (VS) [0.2-0.5 keV], soft (S) [0.5-1.0 keV], 
medium (M) [1.0-2.5 keV] and hard (H) [2.5-10.0 keV]. 
Figure 4 presents the background corrected light curves and 
reveals that the increase around apastron occurs throughout 
the different energy ranges. However, we outline that i) the 
shape of the curves in the H energy band is more symmetric 
around apastron, with an increase that is already present 
at ^ = 0.31 and with the maximum at ^ = 0.54; ii) the 
relative increase is larger for the harder energy bands. The 
observed enhancement indeed corresponds, respectively in 
the S, M and H bands, to about 50, 70 and over 100 per 
cent of the minimal flux level. Such an enhancement of 
the X-ray emission around apastron passage is definitely 
not compatible with a l/Dgep dependence of the fiux, as is 
observed in other colliding wind binaries (e.g. HD 93403, 
Rauw ct al. 2002b) . The corresponding hardness ratios show 
that the emitted X-ray flux is harder around apastron than 
around periastron. Finally, we emphasize that the S band 
accounts for about 50 to 55 per cent of the total detected 
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Figure 4. EPIC MOS and pn light curves (in cnt s~^) of the HD 152248 system in different energy bands. Left panel: 0.2-0.5 keV, 
Middle panel: 0.5-1.0 keV {upper curves) and 1.0-2.5 keV {lower curves). Right panel: 2.5-10.0 keV. Horizontal bars represent the 
effective duration of the exposure while vertical ones axe the 1-a uncertainties on the count rates. 



Table 4. Observed (Obs.) and dereddened (Dered.) X-ray fluxes 
of HD 152248 in the energy band 0.5-10.0 keV corresponding to 
the best fitted 2T-I-PL models (see Table 3). The dereddened 
fluxes were corrected for the ISM absorption only. The corre- 
sponding luminosities were computed assuming a distance mod- 
ulus DM = 11.2. 
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9.32 1032 


0.148 


7.52 10-13 


2.36 10-12 


8.53 1032 


0.312 


7.33 10-13 


2.23 10-12 


8.06 1032 



counts in the [0.5-10.0 koV] band. Spectra of HD 152248 are 
therefore presumably rather soft. 



3.2 Spectra and spectral fits 

As described in Section 2.2, the X-ray spectra of HD 152248 
were consistently extracted for the three instruments 
from each of the six data sets and they were binned to 
reach at least 25 counts per bin2. As shown in Fig. 5, the 
obtained spectra are relatively soft with their maximum 
located between 0.8 and 0.9 keV. At first sight they reveal 
little variability but a general increase of the emission 
level after apastron passage. In order to characterize 
the X-ray spectral properties of HD 152248 and search 
for any modification with time, we investigated several 
spectral models using the xspec v. 11.0 software. In the 
following, we limit our study to energies above 0.5 keV. 
We fixed the interstellar column of absorbing matter to 
a value of nn.isM = 0.311 x 1022 em-2 (obtained from 
Ryter (1996) formula with Ay = 1.4 (Baume, Vazquez 
& Feinstein 1999)). The fitted models are a combination 
of thermal MEKAL models and powerlaw spectra and we 
allow a specific circumstellar column of absorbing matter 
for each component. For each observation, the models were 
adjusted to the individual EPIC MOS and pn spectra. We 
performed simultaneous fits to all three EPIC spectra as 
well. Resulting parameters are reasonably consistcut and 
Table 3 gives the sole results of the EPIC simultaneous 
fitting. 

It became rapidly obvious that a single temperature 

2 Wc also performed the complete analysis using spectra binned 
with at least 10 counts per bin. This docs not alter neither the 
qualitative conclusion, nor the quantitative values of the result- 
ing spectral parameters of the fitted models. However, the ob- 
tained xi were systematically lower than the ones obtained with 
25 counts per bin spectra. 



(IT) model was insufficient to reproduce the observed spec- 
tra {xt > 4). As a second step, we adopted two-component 
models (Fig. 6). We either combined a IT mekal model 
with a powerlaw (PL) or we adjusted two-temperature 
(2T) MEKAL models. Though IT+PL models clearly give 
a better agreement with the data (lower x^)i they fail 
to reproduce the emission lines in EPIC spectra (see e.g. 
the Si XIII multiplet at ~ 1.85 keV in Fig. 6). Restraining 
the fit to the spectral region below 2.5 keV considerably 
improves the quality of the 2T models with respect to the 
IT-fPL models. Actually the power law component helps 
to reproduce the fiux at high energy that is not explained 
by 2T models in some of the spectra that show enough 
signal above 3keV. Adding a PL as a third component to 
the 2T models again increased the quality of the fit for 
those particular spectra (Fig. 7) . The high energy tail could 
also be reproduced by a high temperature (fcTs ~ 4keV) 
MEKAL component. However, except for Obs. 1 {(f) — 0.54), 
a slightly better is reached using 3T models for which 
the three components are soft {kT < IkeV). In this latter 
case however, the hard energy tail might sometimes not 
be adequately reproduced. Though 2T-1-PL models provide 
slightly better fits, it is however tricky to definitely choose 
among all these types of models. If the high energy tail does 
actually correspond to a PL component, this latter is then 
characterized by a photon index P ~ 3.7. This component 
has a much steeper slope than the value P = 1.5 that is 
expected from inverse Compton scattering emission pro- 
duced by a relativistic population of electrons accelerated 
in strong shocks (Chen & White 1991). Table 4 provides 
the observed and dereddened fluxes as deduced from the 
best-fit 2T-I-PL models. 

Though the X-ray flux from HD 152248 deflnitely 
shows a clear variation by about a factor 1.6 between 
its lower and higher value, the spectra reveal very little 
variability at first sight. To detect and to quantify the 
evolution of the fitted spectral parameters in the different 
pointings is not straightforward as different models, but 
also different solutions of the same model, fit the data 
with a similar quality. This latter aspect can be seen from 
Table 3 in which the values of the spectral parameters 
may vary discontinuously from one phase to another. This 
indeterminacy probably results from too large a number of 
free parameters and from the lack of physical constraints 
imposed on the model parameters. We therefore choose to 
limit the fitted interval to the range [0.5 - 2.5 keV] which is 
reasonably well reproduced using a 2T mekal model. This 
approach offers the advantage of avoiding the uncertainty 
about the nature of the hard-energy tail. Though this helps 
to improve the situation, some of the fits still present local 
minima of similar depth. Best fitted values for the column 
of absorbing matter of the first component (nH,i) usually 
tend to be very small, if not equal to zero. Adopting 2T 
models where this parameter is neglected (nH,i = Ocm-2) 
improves slightly but systematically the quality of the 
different fits. It also solves the degeneracy between the local 
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Figure 5. EPIC MOS2 spectra of HD 152248 for the different pointings of the XMM-Newton satelUte. 



Table 3. Spectral parameters of HD 152248 as obtained from three-component models simultaneously fitted to the three EPIC spectra. 
The upper part of this table gives the best-fit 3T models (wabsisM*(wabsi*mekali+wabS2*mekal2+wabS3*mekal3)). The lower part 
reports the best-fit 2T-I-PL models (wabsisM*(w2ibsi*mekali+wabS2*mekal2+wabS3*power)). njj gives the equivalent hydrogen column of 

absorbing matter (in 10^^ cm^^), kT is the temperature of the MEKAL component (in keV) and F is the spectral index of the power-law 

_^ iQ— 14/1 

component. The normalization coefficient of the MEKAL component A'^ (in cm equals ^^^^j J nenB_dV where d is the distance to the 
source (in cm), and rie and nji are the electron and hydrogen number densities (in cm~^), whereas N equals the photon flux at IkeV 
for the PL models. The upper and lower numbers quote the 90 per cent confidence intervals. 
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Two-temperature + powerlaw models 
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minima of the fit, resulting in a coherent fit for the six 

'KMM-Newton pointings. Results of these fits are presented 
in Table 5. The fact that the column of absorbing matter 
of the cold component of the model is ill-constrained could 
be due to the fact that the fitted parameter nn stands for 
an equivalent hydrogen column density of neutral gas. The 
material that is responsible for the local absorption of the 
X-ray flux in the HD 152248 system, however, is ionized due 
to the amount of UV photons emitted by the two O-type 
stars. A neutral absorbing column can then probably not 
properly reproduce the local circumstellar absorption. The 
resulting discrepancies will mainly manifest themselves at 
lower energies where the absorption is the most efficient. 

Figure 8 shows the evolution, with the phase, of 
the best-fit values of the 2T mekal parameters in the 
range 0.5-2.5 keV. The temperatures and normalization 
coefficients of the two components present clear variations. 
These are in fair agreement with the observed modulations 
of the broad band light curves (see Fig. 4). Variations of 
the equivalent column of absorbing matter of the higher 
temperature component are less clear cut. However, there 
might exist two increases at phases (j) = 0.3 and (f> = 0.8 
which roughly correspond to the quadrature phases of the 
HD 152248 system. Though the confidence intervals are 
quite large, this could indicate an enhancement of the local 



Figure 9. Comparison of observed ROSAT PSPC count rates 

(open symbols) and equivalent X.MM- Newton count rates (filled 
symbols) of HD 152248 in the range [0.2-2.4keV]. Phases of 
the ROSAT observations have been computed adopting the 
ephemerides reported in Table 1. 

absorption while the line of sight is almost perpendicular 
to the binary axis and therefore crosses the radial structure 
of the interaction region (sec Sect. 4). 



3.3 Phase-locked variability 

Though the observations clearly show a decrease of the 
emitted flux from Obs. 1-2 to Obs. 6, our XMM-Newton 
campaign alone can however neither establish the phase- 
locked behaviour of the observed variations, nor whether 
this behaviour is stable on longer time-scales. Indeed, 
as seen from Table 2, the XMM-Newton observations of 
HD 152248 extended over about 5 days and therefore did 
not cover more than a single orbital cycle of the binary. 
In order to investigate this question, we retrieved previous 
ROSAT PSPC observations of the NGC6231 cluster 
(Corcoran 1996) and we re-analysed the data adopting the 
ephemerides from Table 1, that are much better constrained 
than the ones Corcoran used at the time. We extracted light 
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Figure 6. EPIC MOSl and pn spectra at </> = 0.54 (upper panels) and </> = 0.15 (lower panels). The spectra were fitted either with 2T 
MEKAL models (lefthand column) or IT+PL models (righthand column) . The model components are drawn with dashed or dashed-dotted 
lines while the solid lines represent the resulting models. Note how the PL component compensates the flux above 3 keV, resulting in a 
fitted model that does not provide a satisfying fit to the emission lines. 



Figure 7. EPIC MOSl and pn spectra at </> = 0.54 (upper panels) and <p = 0.31 (lower panels). The spectra were fitted either with 3T 
MEKAL models (lefthand column) or 2T+PL models (righthand column). The model components are drawn with dashed, dashed-dotted 
or dotted lines while the solid lines represent the resulting models. Note the improvement compared to two-component models. 



curves and spectra in the ROSAT energy band. Due to the 
modest sensitivity and energy coverage of the satellite, the 
PSPC spectra are well reproduced by a single temperature 
component at kT = 0.24 keV and no significant spectral 
variability could be observed in the PSPC data. With 
the help of the xspec software, we then combined the 
2T-I-PL models that best reproduced the XMM-Newton 
observations with the PSPC response file to estimate the 
equivalent PSPC count rates of our XMM-Newton data. 
Figure 9 gives a direct comparison of the equivalent PSPC 
rates with the actual ROSAT observations in the same 
energy range. In view of the difficulty to accurately compare 
the count rates measured from two such different satellites, 
the current agreement between XMM-A^euiton and ROSAT 
observations is quite reasonable and gives a strong support 
to a phase-locked behaviour of the observed variability. It 
also suggests the relative stability, on longer time scales, 
of the dominant phenomenon that produces the observed 
variation in the light curves of HD 152248. Indeed, the 
ROSAT PSPC observations were obtained between March 
1991 and February 1993, i.e. about 10 years prior to our 
XMM-Newton monitoring campaign. 



Table 6. Effective separation and wind velocities at the ram pres- 
sure equilibrium surface on the axis of the system (see text) as a 
function of the HD 152248 orbital phase (p at the time of the six 
yiMM-Newton observations. The last column provides the loga- 
rithm of the total theoretical X-ray luminosity computed for an 
idealized planar collision region (see text). 
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(Sect. 4.1). Sect. 4.2 presents the model used to estimate the 
fluxes that are predicted by the hydrodynamical grids. We 
then report the predictions of the present model (Sect. 4.3) 
and we compare them to the XMM-Newton observations 
(Sect. 4.4). Finally, we conclude by discussing some of the 
obvious limitations of this work (Sect. 4.5). 



4 HYDRODYNAMICAL SIMULATIONS 

As suggested by Pittard & Stevens (1997), comparisons of 
X-ray observations of CWBs with detailed hydrodynamic 
simulations may help to better understand and constrain 
the physics of the phenomenon. Indeed our XMM-Newton 
campaign reveals precious information about the global 
flux emitted by the system that a priori results from 
the X-rays produced in the interaction area as well as 
within the denser layers of the winds of the two O-type 
components. However without further analysis they provide 
little understanding of the local physics of the interaction 
region. Indeed, the 'inversion' of the X-ray data to recover 
the details of the local emission processes as well as the 
various parameters that affect them is certainly a very 
ill-posed problem. In this regard. Computational Fluid 
Dynamics (CFD) has proved to open a new window in the 
physics of the X-ray omission from the winds themselves 
(Feldmeier et al. 1997; Dessart & Owocki 2003) and from 
the wind- wind collision zone (Pittard & Stevens 1997, 2002) . 

This section is organized as follows. First, we give a gen- 
eral description of the code used to solve the hydrodynamic 
problem and we describe the performed computational runs 



4.1 The hydrodynamical code 

In this section, we limit ourselves to a brief overview of the 
hydrodynamical code. A more complete description is given 
by SBP92 and references therein. The VH-1 hydrodynamic 
code aims at the resolution of the Euler differential set of 
equations for inviscid fluids. It is based on the Piece-wise 
ParaboUc Method (Colella & Woodward 1984), a third- 
order accurate time-marching algorithm that is proved to 
be shock-capturing and therefore naturally takes care of 
the shocks. Very briefly, this method is a generalisation of 
the Gudonov scheme that handles each cell interface as 
a shock-tube problem. In this regard, it builds the global 
solution of the equations by solving the local behaviour of 
the fluid (see Anderson 1995, for a general introduction to 
CFD techniques). 

In this work, we adopt the simplified configuration of 
two spherically symmetric winds of constant velocities. We 
therefore both neglect the orbital motion and the wind ac- 
celeration. We assume that the wind material is fully ionized 
(7 = |) and behaves as a perfect gas. Our computations 
account neither for the viscosity, nor for a possible magnetic 
field. These approximations result in an axi-symmetrical 
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Table 5. Spectral parameters of HD 152248 as obtained from 2T mekal models simultaneously fitted to the three EPIC spectra in the 
range 0.5-2. 5keV. Details of the models are: wabsisM*(wabsi*mekali+wabS2*mekal2) in which the value of nH,i has been fixed to zero 
(see text). The same notations as in Table 3 have been used. 



4> 


"H,l 
(cm^^) 


fcTi 
(kcV) 


(10-4 cm-5) 


"H,2 
(10^2 


kT2 

(kcV) 


N2 
(10-4 


xl 


d.o.f 


0.536 
0.659 
0.808 
0.002 
0.148 
0.312 










u.oiUq 

Q 3040.316 
u.outg 289 

u.zyog 284 

0.239 

0.2560:iti 
0.2728:110 


10.12^0797 

io.62li:J| 

ir, CI 10.87 
lU.Ol^o 

9.057|:i? 
8.734|;l|l 
8.127|:f|| 


^.^5*0. 51 

480.53 

Q cqO.63 

410.47 
'^.^^0.34 

n 420.48 

"■""0.53 


r, 71 nO.726 
' ^^0.699 
6940.'''16 

n fiQ40.731 
u.oytQ gg3 

0.6230:148 

0.6210:648 
0.6560:^09 


13 fil 14.24 
ICS. 0^)^2.43 

12.56i3:|i 

9.303j'!,i0 
8.7690:511 
9.183|0g7i 


1.73 
1.51 
1.46 
1.47 
1.45 
1.53 


390 
302 
363 
216 
308 
317 



Figure 8. Best-fit spectral parameters of the 2T mekal models (wabsisM*(wabsi*mekali+wabS2*mekal2)) in the range 0.5-2.5keV 
plotted vs. the phase. Open diamonds refer to the first, low temperature, component parameters while filled diamonds stand for the 
second temperature component (see Table 5). Left: Equivalent column of neutral hydrogen nji. Centre: Temperature kT. Right: 
Normalization coefficient N. 



geometry around the line of centres and allow to reduce the 
hydrodynamical problem to a two-dimensional flow. For 

that reason, all the hydrodynamical runs presented in this 
paper are performed on 2-D grids. We adopt a typical grid 
size of 300 x 600 cells, corresponding to a physical size of 
6 10^^ X lO'^^ cm, and we let the flow evolve with time from 
an initial situation where both winds have not yet collided. 
Wo let the interaction grow while the time elapses and wc 
follow the process for at least 4000 steps, long enough so 
that a steady state may be reached. However, the HD 152248 
CWB harbours a highly radiative collision zone and the 
post-shock flow is dominated by instabilities. Therefore no 
'true' steady state may be reached. Nevertheless the chosen 
evolution time of the simulations is long enough for the 
flow to relax from the initial condition dependencies. The 
time elapsed between two iterations is a free parameter that 
is adjusted at each step according to numerical stability 
considerations, but a typical time interval between two 
steps is about 70 to 120 sec. We performed test runs with 
higher- as well as lower-resolution grids to check that the 
observed behaviour of the flow is not, within the tested 
configurations, resolution-dependent (see however Sect. 
4.5). Energy-loss due to radiation is computed after each 
hydrodynamical iteration and is the only source term 
accounted for in the adopted version of the code. Reflective 
boundary conditions were imposed along the line of star 
centres while we let the gas flow freely out of the grid at 
the other three boundaries. 



As our aim is to compare observations and model pre- 
dictions, the hydrodynamical runs were performed adopting 
a system configuration and wind parameters that closely re- 
produced the configuration of the HD 152248 system at the 
time of the six X-MM-Newton observations. Table 6 reports 
the adopted values for the physical parameters. The separa- 
tions between the stars were computed using the orbital pa- 
rameters and ephemerides tabulated in Table 1. Mass losses 
are deduced from Vink, de Koter & Lamers (2000) formula 
with terminal wind velocities Voo = 2420 km s"^ (Howarth 



et al. 1997) as quoted in Paper I^. HD 152248 being a close 
binary, the winds collide well before reaching their terminal 
velocities. As wind acceleration is neglected, we circumvent 
this problem and adopt the wind velocities reached at the 
position of the ram pressure equilibrium surface on the bi- 
nary axis. This is obtained by solving the on-axis ram pres- 
sure equilibrium equation: 

Mi^^i _ M2V2 ...s 

where di and d2 are the distances from the equilibrium 
surface to either star centre; di + d2 is the separation 
between the centres. Both wind velocities v\ and V2 are 
computed from a standard /J-law with (3 = 0.8. Equation 1 
usually provides three solutions, two of which correspond 
to unstable equilibria. The third one gives the position 
of the stable equilibrium that we then use to derive the 
constant wind velocities reported in Table 6. Though this 
is a rather crude method, it allows to merely bypass the 
problem of non-terminal wind velocities. This provide a 
reasonable approximation to the fact that, as the system is 
rather eccentric, the orbital separation is variable and the 
winds therefore reach quite different velocities at the shock 
surfaces for the different orbital phases that we are studying. 

The hydrodynamical runs performed adopting these 
six sets of parameters finally provide the evolution of the 
values of the gas density, gas pressure and the radial and 
axial velocity components at each grid cell, from which 
other hydrodynamical variables (such as temperature, 
internal energy, entropy...) and radiative properties (see 
Sect. 4.2) can be inferred. Figure 10 gives typical density 

3 In Paper I, we reported a slightly erroneous value for the termi- 
nal velocity of the secondary component. The present value (see 
Table 1) is therefore corrected accordingly. Fortunately, this does 
not alter the mass loss rate computation. Our previous conclu- 
sions remains therefore unchanged except that, using the present 
corrected values, the on axis wind momentum ratio is now slightly 
in favour of the secondary star. 
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Figure 10. Density (left) and temperature (riglit) maps as computed from our liydrodynamical simulations. The adopted parameters 
(see text) reproduce the HD 152248 configuration at phase </> = 0.15. 



and temperature maps provided by the simulations of the 
HD 152248 CWB system. 



Figure 11. Cross-section at z = 2 10^^ cm in the density (dashed 
line) and temperature (solid line) maps presented in Fig. 10. Both 
variables are plotted using logarithmic scales. 



4.2 X-ray emission 

Prom the results of the simulations described in the previous 
section, we extracted maps of the hydro dynamical configu- 
ration (i.e. density, pressure, temperature, r- and z-velocity 
component maps ...) once every 10 time steps. Depending 
on the initial configuration, it takes about 2100 to 2600 
iterations for the system to relax from the initial conditions, 
which corresponds to an elapsed time ranging from about 
1.5 to 2.0 10* sec. Once the dynamical flow has relax;ed, we 
obtained a snapshot approximately every 700 to 750 sec. 
From the succession of these snapshots, we thus follow 
the evolution of the hydrodynamical variables with time, 
and therefore of the wind-wind interaction phenomenon, 
in the frozen configuration of the HD 152248 system that 
corresponds to one of the six XMM-Newton observations. 
By themselves, the hydrodynamical simulations are very 
instructive and it is worth to monitor the flow structure of 
the collision area. However, to compare these simulations 
with observational results, we need to deduce observable 
variables from the hydrodynamical grids. For this purpose 
we then solve the radiative problem corresponding to the 
different hydrodynamical conflgurations. 

The X-ray emission from the system was evaluated by 
summing up the emissivity of each cell, accounting for a 
proper extinction (see details in SBP92). The computation 
of the column of absorbing matter along the line of sight 
required a three-dimensional geometry as we have to 
account both for the inclination of the system and for 
the variable line of sight of the observer. A 3-D grid was 
thus generated by a rotation of the 2-D hydrodynamical 
grid around the line of centres. The adopted angular 
resolution of the constructed 3-D cells is 2°. The chemical 
composition of both winds was supposed to be solar. For 
each cell, we computed the emissivity as a function of the 
mean cell temperature. The contribution of grid cells with 
a temperature lower than 10® K was neglected. For the 
emitting cells, the corresponding column of circumstellar 
absorbing matter towards the observer was computed and 
we used proper absorbing coefficients to derive the intrinsic 
attenuated spectra. Cells with temperature higher than 
10^ K were considered not to give a significant contribution 
to the absorption. Eclipses were accounted for assuming 
stellar radii of 15 R©. Interstellar absorption was taken 
into account assuming a column of absorbing matter of 
?^H,ISM = 0.311 X 10^^ cm~^. We finally obtained predicted 
X-ray fluxes and spectra that are suitable for a direct 
comparison with the X-MM-Newton observations. 



4.3 Predictions of the simulations 

As can be seen from Fig. 10, the collision zone turns out 
to be highly unstable. The wind-wind interaction region 
appears to be extremely thin and to develop filaments 
that form smoke-like structures while time elapses. Similar 
patterns arc found for the six configurations in which we 
performed the simulations. The thinness of the collision 
zone is expected due to the highly radiative behaviour of the 
collision (x = 0.02 — 0.15). The interaction region appears 
to be so narrow in our simulations that its inner structure 
could not be resolved with our mesh size. Increasing the 
resolution by a factor two does not solve the problem. This 
failure to resolve the collision area constitutes an obvi- 
ous limitation of our work that will be discussed in Sect. 4.5. 

The filaments that form the interaction region slowly 
evolve to produce thin but dense shells that axe separated 
by low density regions (Fig. 11). In addition to this general 
pattern, some pockets of high temperature gas seem to 
undergo an adiabatic expansion with typical cooling time 
scale larger than a couple of weeks. It is probable that 
some of these pockets might survive till far from the axis 
but, due to their very low density, it is unlikely that they 
provide a significant contribution to the total emission from 
the wind-wind interaction. 

To investigate the collision zone luminosity, we com- 
puted the expected emission for each of the extracted 
snapshots. One of the surprising results of these simulations 
is that the luminosity from the wind collision can vary 
by a factor of a few on time scales of a couple of hours 
(see Fig. 12). As the time resolution of our simulations 
is about 750 sec, we also performed a run in a particular 
configuration, extracting snapshots after each hydrody- 
namical iteration. This provides a light curve with a time 
step of about 75 sec. The latter almost exactly matches the 
previous light curve that has a resolution ten times lower. 
We can therefore reasonably assume that all the temporal 
information concerning the variability of the luminosity 
are already contained in the light curves presented in Fig. 
12. To estimate the expected luminosity at a given phase 
of HD 152248, we therefore computed the average and 
dispersion of the luminosity of the selected snapshots. This 
is presented in Fig. 13 where the vertical bars should be 
considered as the model envelope rather than as 'classical' 
error bars. The upper panel gives the variation of the 
intrinsic emission produced in the collision of the winds. 
The latter displays a variation of about 40 per cent between 
periastron and apastron. However accounting for the 
eclipses and the absorption by the circum- and interstellar 
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Figure 12. Evolution witli the time of the logarithm of the ab- 
sorbed X-ray luminosity emerging from the interaction region in 
the range 0.5-10.0 keV, as predicted by the models. Different lines 
refer to different configurations of the system. Solid line: <j) = 0.54. 
Dashed Une: </> = 0.00. Dotted Une: </> = 0.15. 



Figure 13. Simulated X-ray luminosity from the collision zone 
in the range 0.5-10.0 keV as predicted by the models at the dif- 
ferent phases of HD 152248. Upper panel: Open circles: Intrinsic 
emission of the collision zone. Filled triangles: theoretical X-ray 
luminosities from the collision of two identical winds (see Table 6) . 
Lower panel : Emerging luminosity accounting for i) eclipses, ii) 
local absorption, iii) ISM absorption. The vertical bars represent 
the 1-(T dispersion on the emerging fluxes. 



medium drastically changes the relative amplitude of the 
variation. Prom the lower panel of Fig. 13, it is obvious 
that the performed hydrodynamical simulations predict a 
phase-locked variation of the outgoing flux by about one 
order of magnitude. The maximum of the flux is reached 
around apastron, while the emitted flux is minimum at 
periastron. 

In order to check the quality of our model, we also com- 
puted the theoretical luminosity emitted by the collision of 
two winds of equal strength that would form an idealized 
planar interaction region. In this simple model, we further 
assumed that the kinetical energy of the flow component ra- 
dial to the shock surface is fully converted into X-rays. For 
each wind, this is given by (Luo, McCray & Mac Low 1990): 



7-th _ Lyi^i 

J-i-x i — — :: — 



where L„ 



12 



(2) 



: — is the equivalent wind luminosity of 
star i, characterized by a mass-loss rate Mi and a pre-shock 
velocity Vi. Typical values for the kinetic energy are of 
the order of 10^® — 10^^ erg for single O-stars with 
terminal wind velocities. Values reported in Table 6 for the 
HD 152248 system result from the sum of the X-ray lumi- 
nosities from the two shocked winds : Lx = -'^x,! + ix,2- 
These values are in excellent agreement with the averaged 
intrinsic luminosities predicted by the hydrodynamical 
simulations (see Fig. 13) and lend further support to the 
consistency of our numerical results. 



4.4 Comparison with the XMM-Newton 
observations and discussion 

The hydrodynamical simulations of wind-wind collisions 
performed in configurations as close as possible to the 
ones of the HD 152248 system at the different phases of 
interest predict an increase of the emitted fluxes from the 
interaction region by about one order of magnitude around 
apastron. Though the observed EPIC light curves indeed 
reveal an enhancement of the emission at this particular 
phase, the increase is limited to a factor of about 1.6 
and further shows a clear asymmetry between the phase 



Figure 14. Comparison between the computed dereddened lu- 
minosities as resulting from our simulation runs and the obser- 
vational dereddened luminosities. Filled squares: dereddened lu- 
minosities of the interaction region as predicted by the model. 

The vertical bars have the same meaning as in Fig. 13. Filled cir- 
cles: total predicted luminosities, i.e. resulting from the predicted 
emission of the interaction region plus the expected contribution 
from the two components of HD 152248 (see text). Open circles: 
observational dereddened luminosities as deduced from the fitting 
of the XMM EPIC spectra (see Table 4). 



interval [0.5-1.0] and [0.0-0.5]. Due to the 2-D nature of the 
simulations, the predicted light curve is naturally symmet- 
ric around apastron and, in this regard, does not render 
the observations. This might suggest that the observed 
asymmetry could partly find an explanation in the orbital 
motion and the resulting deflection of the colliding wind 
zone. 

Assuming a distance modulus DM = 11.2 (Raboud et 

al. 1997) as the one we used to infer the bolometric lumi- 
nosities and subsequent mass-loss rates reported in Paper I 
and Table 1, the predicted model luniiuosity clearly under- 
estimates the 'true' luminosity of the system (see Table 4). 
However the luminosities computed from the hydrodynam- 
ical simulations only account for the emission coming from 
the wind-wind collision itself. Until now, the intrinsic X-ray 
emission of the two O stars has indeed been neglected. From 
Paper I, we know that log(L]iJ"™/L0) and log(L^o'i/-^0) 
equal 5.61 and 5.63 respectively for the primary and sec- 
ondary component. Using the relation of Berghofer et al. 
(1997) for stars with Lboi > 10^** erg s"^: 



log (L> 



1.13 log (Lboi) - 11.89 



(3) 



we find that the intrinsic X-ray luminosities of the com- 
ponents of HD 152248 should be about log(LP"'") = 32.39 
and log(Lf'' ) = 32.42 (erg s~^) respectively. Adding the 
dereddened flux originating from the interaction region 
therefore provides a total X-ray luminosity suitable for 
a comparison with the observed values. We outline that 
the Berghofer et al. relation is based on observed X-ray 
emissions from O-stars and so includes the effect of in- 
trinsic wind absorption. The resulting luminosities of the 
components of HD 152248 are however not corrected for 
the additional absorption associated with the colliding 
wind region nor for the eclipses that occur in the system. 
Figure 14 presents a direct comparison of the dereddened 
observational X-ray luminosities with the luminosities 
predicted by the simulations. The latter are obtained by 
summing up the emission from the interaction region with 
the expected intrinsic contributions from both stars of the 
system. Though there is still a difference of about a factor 
two between the observations and the model for some of the 
pointings, the agreement is now much better. The intrinsic 
dispersion of the Berghofer et al. relation and the resulting 
uncertainty on the intrinsic Lx can, by itself, be responsible 
for such a discrepancy. The amplitude of the variation 
between apastron and periastron in the model and in the 
observed light curve differs by about 50 per cent. Given the 
approximations inherent to our simulations, the agreement 
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Figure 15. Comparison between the computed spectra as result- 
ing from our simulation runs and the observational best-fit model 
for the configuration of the HD 152248 system at <^ = 0.31. Plain 
thin line: intrinsic predicted spectrum. Dashed line: absorbed pre- 
dicted spectrum. Vertical bars stand for the model envelope due to 
fluctuations in the predicted spectrum. Thick line: best-fit 2T-I-PL 
model at (/) = 0.31 (see Table 3). 



can be considered as quite good. 

Though the quality of our model is hmitcd (sec Sect. 
4.5), these results suggest that the observed phase-locked 
behaviour of the X-ray flux of HD 152248 may indeed be 
reasonably explained by a variation of the X-ray emission 
produced in the interaction region. While the intrinsic 
colliding wind flux increases at apastron duo to larger pre- 
shock wind velocities, the major effect is to be attributed 
to the variation of the absorption properties. The predicted 
modulation, if diluted by the intrinsic emission of the two 
components of the system, indeed reasonably reproduced 
the observed amplitude of the fluxes deduced from the 
XMM EPIC spectra. 

The energy distribution of the predicted X-ray emission 
presents its maximum around 3kcV. Figure 15 provides 
a direct comparison with the best-fit model 2T-I-PL that 
reproduces the observational spectral properties of the 
HD 152248 system. Similar patterns are found in the six 
configurations. It is seen that the predicted spectra are 
much harder than the observed ones and, on average, 
underestimate the flux at low-energy (below 2keV) while 
the flux is overestimated at higher energies. However, as 
previously emphasized, the predicted spectra only account 
for the wind-wind collision emission and the X-ray fluxes 
from the two stars are neglected. The X-ray emission 
from individual stars is expected to be relatively soft, 
with its maximum peak being reached around 1 keV. This 
can therefore account for the observed discrepancy at 
low energy. On the other hand, the disagreement of the 
high-energy tail might come from an excessive strength of 
the shocks in the model, possibly resulting from too large 
pre-shock wind velocities being assumed. Wind acceleration 
might indeed be slowed down by speciflc effects such as 
radiative iuhibition (Stevens & Pollock 1994). However, the 
boundaries of the spectra predicted by the hydrodynamical 
models still overlap the observational best-fit model at 
high-energy. 

Finally, a moderately strong 6.7 keV Fe line is clearly 
predicted by the simulations while its presence is not so 
obvious in the data presented until now. The observational 
model plotted in Fig. 15 is a 2T+PL model in which the 
PL component accounts for the high energy tail and could 
intrinsically not reproduce the Fcxxv line at 6.7 keV. De- 
termining whether this line is present or not in the spectrum 
of HD 152248 might help to clarify the nature of the hard 
energy tail. For this purpose, wo combined the EPIC data 
of the six X.MM-Newton pointings to increase the signal- 
to-noise ratio at high energy. We built merged event lists 
and created combined MOS and pn arf using the ftools 



Figure 16. Merged EPIC pn spectrum of the HD 152248 binary. 
The spectrum is plotted using 17 counts per bin for the sake of 
clarity. The inset presents the result of a PL + Gaussian model 
adjusted in the vicinity of the 6.7keVFe line (see text). 



command addarf. The merged pn spectrum of HD 152248 
is presented in Fig. 16 and corresponds to a total exposure 
time around 130 kscc. We then adjusted an empirical model 
(power+gauss) in the range [4.0-10.0 keV]. This model is 
made of a PL component, that reproduces the continuum, 
and of a Gaussian component (-^7|= — exp (— ( )^) ) 
that fits the Fe line (see Fig. 16). Best- fit values for 
the Gaussian parameters give -Eg = 6.70 ± 0.14keV, 
CTG = 0.24 ± 0.14 keV and /sT = 4.3 ± 2.6 10"'^ cnt s-^cm-^ 
yielding an equivalent width of about 2.3 keV. From the 
merged data, we can infer that the 6.7kcV Fe lino is 
most probably present in the HD 152248 spectrum and 
that, therefore, at least part of the high energy tail has 
a thermal origin. Finally, we mention that De Becker et 
al. (2004b) recently noted that non-thermal phenomena 
could also contribute to the Fexxv emission at 6.7 keV. 
The above diagnostic seems however to remain quite robust 
as these authors mentioned that the relative contribution 
of the non-thermal and thermal processes to the emission 
is directly proportional to the fraction of the free electron 
population that reaches rclativistic energies. It is therefore 
to be expected that a non-thermal contribution to the 
Fexxv line remains marginal. 



4.5 Limitations of the model 

A first obvious limitation of these simulations is the use 
of constant wind velocities in a geometry for which we 
definitively know that the collision occurs in the accelera- 
tion region of the winds. Adopting the velocities reached 
at the ram pressure equilibrium surface offers nevertheless 
the advantage of a quite straightforward solution. This 
approach probably provides a good approximation close to 
the line of centres. However, our simulations reveal that the 
structure of the interaction extends far from the binary axis 
and the hypothesis of constant wind velocity is definitely 
too simplistic there. One of the ways to account for the 
line-driven acceleration would be to include an extra source 
term that would express a radiative force on the wind 
material. Due to the wind acceleration, the pre-shock wind 
velocities away from the axis would definitely be larger and 
the shocks stronger. Under the influence of an increased 
ram pressure, aud of its component normal to the shock 
surface, the interaction region could be more confined along 
the radial direction. However such an effect is difficult to 
estimate and, furthermore, would probably be disrupted by 
the orbital motion in the system. 

Moreover, as seen from Fig. 10, our grid is too coarse to 
resolve the inner structure of the interaction region. Though 
a thin interaction is expected duo to the importance of 
the radiative cooling, a linear increase of the resolution 
would however not have helped to solve the problem. 
Indeed we are facing the need to simultaneously resolve 
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small scale structures around and inside the interaction 
region, but also to follow the evolution of large scale 
structures such as the shells. However computation time 
prevents us from increasing both the resolution of the 
mesh and simultaneously maintaining a good coverage 
of the collision far from the axis. From the density and 
temperature conditions of the interaction region inferred 
from our simulations, we estimate the typical cooling length 
of the shock-heated plasma to be around 10^ — 10* cm. 
On axis however, extreme values of a few 10® cm can be 
found. Our grid resolution along the axis is 2 10^" cm and 
is definitely not designed to resolve the inner structure of 
the collision zone. Increasing our resolution by a factor of 
10 , is probably required to overcome that problem. This 
is unfortunately not realistic with the present hardware 
resources. In such a situation, the need to adopt non-linear 
metrics that would increase the resolution towards the 
central axis of the systems is obvious. An adaptative grid 
that automatically clusters grid points in the region of high 
flow-field gradients (see e.g. Anderson 1995) certainly offers 
another powerful way to overcome this problem. Due to 
the instability of the region and the complex structure that 
results from it, the mesh should also ideally be readjusted 
after each iteration. Therefore, to adequately perform 
hydrodynamical simulations of wind- wind collisions in close 
binary systems, an adaptative mesh algorithm, able to 
identify the regions where an increased resolution is needed, 
can provide an adequate description using reasonable time 
and resources. The implementation of such techniques is 
not trivial, but remains one of the most promising way to 
drastically increase the quality of numerical simulations 
of close highly radiative CWBs. Nevertheless, as seen 
from Fig. 13 and Table 6, the macroscopic properties, and 
more specifically the intrinsic emission, of the interaction 
region in our simulation are in fair agreement with what is 
expected from purely theoretical considerations. This lends 
good support to our results, even if it is clear that a linear 
grid as the one we used is certainly not the most optimized 
approach to the current problem. 

Other limitations are related to the 2-D nature of the 
simulations that inherently neglects the orbital motion of 
the binary. In a system such as HD 152248, the orbital 
velocities of a few hundred km s~^ compete with the 
pre-shock wind velocities that are only 2 to 4 times larger. 
Under such circumstances, the hypothesis of isotropic winds 
might also break down. Furthermore, in our simulations, 
the shell structure has radial progression velocity of only 
100 km s~^, much slower than the orbital motion. It is 
therefore very likely that the orbital motion might have a 
considerable influence on the geometry of the interaction 
region. One of the most obvious effects is the deflection 
of the interaction region. Due to the competition between 
the orbital motion and the radial velocities of the shell 
structure, we expect the interaction zone to be rather 
tilted. Another problem inherent to the 2-D nature of the 
simulation is that we let the flow evolve in a frozen orbital 
configuration. Therefore the X-ray emission obtained is 
definitely independent of the history of the system. The 
winds, escape flow and orbital motion of the HD 152248 
system seem to evolve on similar time-scales. In the real 
world, one can therefore reasonably expect that these 



phenomena compete or combine and that the resulting 
X-ray emission is correlated to the history of the different 
phenomena. We also note that tidal distortions of the stars 
in an eccentric binary may induce an asymmetric behaviour 
of their interaction with respect to periastron (Moreno & 
Koenigsberger 1999). 

To conclude this section, the physics of the model 
itself can still be improved. Notably by including in the 
sinnilations other physical phenomena such as radiative 
braking (Gayley, Owocki & Cranmer 1997) and/or inhibi- 
tion (Stevens & Pollock 1994). In the particular case of the 
HD 152248 binary, sudden radiative braking is unlikely to 
play a major role as both stars of the system are of a broadly 
similar nature and therefore display similar luminosities and 
have winds of similar mass-loss rates and velocities. These 
winds will thus have similar radiative driving efficiencies 
(that is similar CAK type wind constants k and a). All 
these factors work against sudden radiative braking being 
important in this system and we can safely ignore it as 
a major effect for the studied object. A potentially more 
significant problem is the inhibition effect. Since the stars 
arc quite close, the radiation of both stars will play an im- 
portant role and could alter the geometry of the wind-wind 
collision, both in terms of slowing and obliquely deflecting 
the colliding wind stream. Finally, the influence of the X- 
rays emitted by the collision on the winds of the stars and 
tidal distortions of the star surfaces that probably affect the 
structure of the winds are further effects that should also be 
addressed in CWB simulations. Properly accounting for all 
the above mentioned physical phenomena within the code is 
however far beyond our present purposes and remains part 
of long term developments of CWB hydrodynamical models. 



5 SHORT TERM VARIABILITY 

As suggested by the hydrodynamical simulations, the 
X-ray light curves of HD 152248 might display evidence 
of short-term variability, i.e. variations on characteristic 
time scales ranging from a couple of minutes to a couple 
of hours and that are therefore much shorter than the 
orbital period. This section summarizes the results of the 
different methods implemented to search for short term 
variability within the collected X-MM-Newton data. For 
each pointing, we applied the following variability tests: 
i) we fitted polynomials of degree n, with n = 0, 1 or 2, 
and we estimated both the quality Q of the fit and the 
merit of including an additional term; ii) we applied 
the Kolmogorov-Smirnov test to check the uniformity of 
the distribution function; iii) we used a corrected version 
of the probability of variability test suggested by Preibisch 

6 Zinnecker (2002, sec Appendix A); iv) finally, we also 
performed a Detrended Fluctuatton Analysis (DFA), a 
method that aims at searching for long-range correlation 
within time series. The latter three methods were applied 
on the series formed by photon arrival times. Though it 
was therefore not possible to account for a background 
correction, we restrain our analysis to PI events larger than 
500 to minimize the impact of background events. Indeed, 
in the [500-10 000] PI energy band, the background level 
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Figure 17. Background corrected light curve (in cnt s ^) of the 
HD 152248 system in the range PI e [500 - 10 000] as presented 
in Fig. 3. Detected local trends within the different sets of obser- 
vation (see Table 7) have been overplotted (plain lines). Dashed 
lines were computed assuming a l-cr difference from the best-fit 
slope values. The fitted local trends of Obs. 2 and 6 (at respec- 
tively rp = 0.66 and 0.31) arc not significantly different from a 
constant level and are therefore not plotted on this figure. 



Table 7. Best fit values of the slope whenever a lineaj: model 
is the most appropriate to describe the light curve. The adopted 
bin size is 1000 sec and the energy ranges approximately from 0.5 
to lOkeV (Pie [500- 10 000]). All the values are expressed in lO"'' 

cnts~^. 



EPIC 


Obs. 1 


Obs. 3 


Obs. 4 


Obs. 5 


MOSl 


9.6 ±2.4 


-5.4 ±2.1 


-4.5 ± 2.1 


4.7 ±2.2 


MOS2 


5.3 ±2.4 


-8.4 ±2.1 


-2.5 ± 2.2 


7.6 ± 2.2 


pn 


11.0 ±4.0 


-19.5 ±3.6 


1.2 ±3.8 


11.6 ±3.7 



only accounts for about 2% of the detected counts and 
should therefore not be a critical issue for these tests. 



5.1 Polynomial fits 

For each pointing, we first extracted background corrected 
light curves from the filtered event list using time bin sizes 
ranging from 10 sec to 5000 sec in order to investigate the 
difi^erent time scales. We then fitted polynomials of degree 
n = 0, 1 and 2 to the light curves. In doing so, we only 
consider here the detection of trends within a single 30 ksec 
observation. The following criteria were used to estimate 
the adequation of the fitted model. First, we computed the 
merit of adding a supplementary term to a polynomial 
of degree n. This is given by : 



■ Xn+1 



xl+J{N-n-2) 



(4) 



where N is the number of bins in the light curve (see e.g. 

Bcvington 1969). As a ratio of two distributions, the 
statistic F-^ follows the Snedecor F-distribution with v\ = 1 
and V2 = A*' — n — 2. We then required that F-^ exceeds a 
threshold value F-^ > /, corresponding to the probability 
Pf(/; i^i, 1^2) = 0.95. Second, the model should reasonably 
fit the light curve, i.e. the quality Q of the fit - which is 
given by: 

Q(^,^) = l-^^(^,^) (5) 

where P is the incomplete gamma function and v, the 
number of degrees of freedom - should yield a reasonable 
value (see e.g. Press et al. 1989). Table 7 reports the 
values of the trends whenever a linear model was the most 
appropriate to describe the data. Though these values seem 
to be quite small in terms of count rate, they might account 
for a 10 to 20 per cent variation over a 30 ksec exposure 
and are therefore significant. The inclusion of an additional 
square term turns out to be unnecessary in most of the 



cases and only marginally improves the quality of the fit 
in the very few cases where an additional term could have 
been relevant. These constraints on the values of the first 
order derivative of the flux might be useful in testing more 
evolved models where, for example, the orbital evolution of 
the system with time is accounted for. 

We note that the slopes observed for the two MOS 
instruments are in reasonable agreement, whereas the pn 
best fit slope is usually steeper. This might be explained 
by the different response of the MOS and pn and, indeed, 
the slopes observed in the S and M energy bands are often 
significantly different. This is further in agreement with 
the fact that the global variation of the fiux goes along 
with a modification of the spectral properties of HD 152248 
(see Sect 3.2). Figure 17 displays the background corrected 
lightcurve of Fig. 3 on which the detected local trends have 
been represented. It is clear from this figure that most of 
the detected trends are indeed in line with expectations 
from the longer timescale orbital variation, but for Obs. 
5 at = 0.15. Indeed the steep trend detected in Obs. 5 
count rates seems to initiate a rising of the flux towards 
apastron, as it is seen in the model predictions. However, 
Obs. 6 at <^ = 0.31 does clearly not follow that trend and 
shows a further decrease of the flux. The exact reason of the 
limitation of the fiux increase towards apastron remains at 
this stage largely unexplained. However the simulations of 
the previous section indicate that part of the answer might 
be linked to the orbital motion of the system. 



5.2 Kolmogorov-Smirnov 

To check for constancy, the Kolmogorov-Smirnov test was 
applied on the series formed by the time difference between 
two photon arrival. We accounted for the Bad Time 
Intervals (BTIs) within a single pointing by subtracting the 
respective BTI durations to the arrival time of the photons 
that were chronologically following the BTIs. This is equiva- 
lent to assuming that the detector is switched off during the 
BTIs. Such a procedure might actually favour the rejection 
of the null hypothesis of a uniform distribution if and 
only if the underlying distribution is not uniform. In this 
sense, we might therefore consider that the adopted elimina- 
tion of the BTIs does not bias the Kolmogorov-Smirnov test. 

The Kolmogorov-Smirnov test reveals significant 
variability (> 95 per cent) for those pointings that display 
trends, as shown in the previous paragraph (Sect. 5.1). 
Inspection of the time series further allows to assert that 
the Kolmogorov-Smirnov test actually detects the local 
trends within the relevant pointings, in agreement with the 
results of the polynomial fits (Sect. 5.1). 



5.3 Probability of Variability 

We applied a modified version of the Prcibisch & Zinnecker 
(2002) probability of variability test (see Appendix A) 
to the HD 152248 and background event lists for which 
the PI events are in the range [500-10 000[. One of the 
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Figure 18. Integrated time scries y{i) = ('^'j ~ built 

from the EPIC pn observation at (j> = 0.81 in two energy ranges. 
Rl; Pie [500 - 10 000] (solid line) and R3: Pie [2000 - 10000] 
(dotted line). Left and bottom scales refer to Rl while right and 
top scales refer to R3. 



Figure 19. Evolution of F{k) with the window size k plotted 
on a log — log scale diagram. The data set is built from EPIC 
pn events at </> = 0.81. Different symbols refer to different energy 
ranges. Notations Rl to R3 have the same meaning as in Fig. 18 
and in Table 8. 



advantages of this method is to directly account for the 
poissonian distribution of low count rate observations 
and it is therefore particularly designed to investigate any 
background variability. As the values of the adopted binning 

depends upon the expected number of counts < A*' >, wc 
choose < A'' > to range from two to a few thousands in 
order to investigate different time scales. Indeed a larger 
value for < N >, and licucc for the time bin size, favours 
the detection of long range variation. On the other hand, 
small values of < A'^ > will tend to enhance small scale 
variations that could be diluted or averaged out in the first 
case. The choice for < iV > could thus enhance or hide 
variability within a particular characteristic time scale. We 
adopt a confidence level of pov > 0.99 for the rejection of 
the mill hypothesis of a uniform distribution of the counts 
throughout the exposure. 

Concerning the HD 152248 binary system, clear cut 
positive results occur for Obs. 3 (i^ = 0.81) for the EPIC 
MOS2 and pn instruments while variability in the MOSl 
is only marginally detected (0.90 < pov < 0.99). The 
variability seems to occur on time scales above 1000 sec 
and might be related to the detection of a 'steep' trend 
in the related data. The background of Obs. 4 (0 = 0.00) 
is also variable, especially for the EPIC pn instrument 
though the corresponding HD 152248 event list for Obs. 4 
gives a negative result under the pov test. Finally marginal 
variability (0.90 < pov < 0.99) is also detected for Obs. 1 
{4> = 0.54) and might again be related to the detected trend. 

In conclusion, the three methods described in Sections 
5.1 to 5.3 give fully consistent results and prove that 
local trends, presumably linked to the orbital variability, 
are clearly present within four out of the six pointings. 
However, shorter time scale variability such as those seen 
in the results of the models described in Section 4 has 
not been found. We might wonder if the local trends do 
not alter the properties of the time series and therefore 
obscure the possible intrinsic variability of the wind-wind 
collision. One of the recently developed methods that allows 
to characterise a signal intrinsic variability regardless of 
possible local trends is the Detrended Fluctuation Analysis. 
This will be the focus of the next section. 



5.4 Detrended Fluctuation Analysis 

The Detrended Fluctuation Analysis (DFA) is a method that 
aims at searching for long-range correlations in noisy signals 
(Peng et al. 1994, 1995) and is based on the variance analysis 
of the fluctuations around the local trend. Assuming a self- 
similar sequence {y{i)} of A'^ elements, we divide the series 
into A^w non-overlapping windows of A: = N/N^ elements. 
We then compute the mean square fluctuation around the 
local trend yi in the considered window : 

(l + l)k 
i=lk+l 

where yi is the least-square linear fit function in the window 
I {I = 0, 1,2, A^w - 1). We then average Fi{l) over the 
A'w windows of size k. Repeating the operation for different 
k values yields a power-law relation between the root mean 
square fluctuation function F and the window size k (see 
e.g. Kantelhardt et al. 2001): 

m = V^n^ - (7) 

where a is the Hurst exponent, also called the long-range 
correlation exponent, which is typical of the series be- 
haviour. For a purely random walk (e.g. Brownian motion), 
a = |. For a < I, the signal harbours anti-persistent 
correlations while for a > i, persistent correlations are 
present (see e.g. Feder 1989). This method has proved its 
worth in a wide variety of domains such as physiology 
(Ashkenazy et al. 2001), biophysics (Stanley et al. 1993) 
and econophysics (Vandewalle & Ausloos 1997). Stanley 
et al. (1993) also compared its performances to other 
more widespread techniques (correlation function, power 
spectrum, rescaled range (R/S) analysis) and report that 
the performances of the DFA method merely rest on its 
increased ability to reduce the noise in the determination of 
the correlation exponent a and to the fact that the method 
is mostly insensitive to both local and global trends. Some 
authors (e.g. Hu et al. 2001) suggest the use of polynomials 
of degree n (DFA-n analysis), with n larger than one, 
instead of the linear local trend yi. We however already 
showed (Sect. 5.1) that linear models are the most successful 
in describing the trend within a single pointing. We there- 
fore limit our DFA analysis to the so-called DFA-1 approach. 

As the extension of the method over several decades of 
the window size k is crucial to properly estimate the cor- 
relation exponent, we construct our time series using the 
time intervals Ati that separate the arrivals of two succes- 
sive photons in the adopted region for HD 152248 (see Sect. 
2.2). The constructed time series are then mapped onto a 
self-similar process by integration: 

i 

j/(i) = 5I(Ai.-At) (8) 

i=i 

where At is the average of the time intervals At over 
the whole time series. Figure 18 shows an example of 
such a series corresponding to the EPIC pn observation 
of HD 152248 at </> = 0.81. For a window size ranging 
approximately from fc = 4 to = A^/4, we then compute 
the average F{k) values. As can be seen from Fig. 19, this 
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Table 8. Values of the Hurst exponent (or long-range correla- 
tion exponent) a as obtained from a least-square linear fit on a 
log — log scale for the three EPIC instruments and for different 
energy ranges. Rl: Pie [500 - 10 000]; R2: Pie [1000 - 10 000]; 
R3: Pie [2000 - 10 000]. 



Obs. # 


PI 


MOSl 


MOS2 


pn 


1 


Rl 


0.540±0.006 


0.540±0.004 


0.512±0.004 




iXZ 


U.OozztU.UUy 


\J.D^( ztU.UUD 


n c;i Q_i_A nrxfi 




R3 


0.555ib0.016 


0.535ib0.016 


0.537ib0.008 


2 


Rl 


0.496±0.008 


0.495±0.009 


0.506±0.006 








n 0-1-0 n07 


n ^4-0 noR 




R3 


0.470±0.021 


0.681±0.056 


0.674±0.019 


3 


Rl 


0.502±0.009 


0.500±0.007 


0.514±0.007 




R2 


0.479±0.014 


0.518it0.012 


0.503±0.006 




R3 


0.563±0.021 


0.589±0.025 


0.568±0.014 


4 


Rl 


0.521±0.009 


0.496±0.009 


0.546±0.005 




R2 


0.496±0.011 


0.503it0.012 


0.514±0.013 




R3 


0.549±0.021 


0.583±0.025 


0.517±0.022 


5 


Rl 


0.553±0.011 


0.521±0.009 


0.502±0.004 




R2 


0.497±0.011 


0.495±0.011 


0.461±0.007 




R3 


0.535±0.021 


0.687±0.027 


0.496±0.017 


6 


Rl 


0.548±0.007 


0.477±0.008 


0.500±0.005 




R2 


0.521±0.013 


0.475±0.010 


0.508±0.007 




R3 


0.547±0.019 


0.494±0.022 


0.514±0.015 



yields a power-law relation F{k) ~ k". We then adjust a 
linear relation on a log — log scale to estimate the value of 
the power-law parameter a. Figure 19 presents a typical 
evolution of F{k) with the number k of elements in a 
window and Table 8 gives the values of the Hurst exponent 
a as obtained from a least-square linear fit. We performed 
a similar analysis in different energy bands. While in the 
range PI e [500 — 10 000], a takes values around 0.5, it 
appears that the self-similarity parameter a tends to reach 
somewhat higher values when restricting the data set to 
higher energies. This suggests that persistent long-range 
correlation might be present in the higher energy tail of the 
flux. As indicated by our simulations, this hard tail above 
2 keV might find its origin in the interaction region. 

Sreenivasan &: Meneveau (1986) showed that several 

facets of turbulence could be related to fractals. In partic- 
ular, Meneveau & Sreenivasan (1987) reported that energy 
dissipation in turbulent processes displays long-range 
correlations with a < 0.5. Their best model yields a « 0.45. 
The apparent increase of the self-similarity coefficient with 
the energy in our XMM-Wewton time series tends to yield 
a values slightly larger than 0.5 and might therefore not 
be directly related to turbulence occurring within the 
wind-wind interaction region. Though the results of the 
DFA method presented in this paper are not clear cut, they 
however constitute an encouragement to apply this kind 
of technique to a wide variety of astrophysical time series 
analyses. 



6 SUMMARY AND CONCLUSIONS 

In this paper, we presented the results of an XMM-Newton 
campaign on the massive CWB binary HD 152248 that 
was performed as part of the GT observation time of the 
OM-consortium. The campaign was split into six separate 
pointings, for a total duration of 180 ksec and the phase res- 
olution of the observations is probably one of the best ever 
performed on an 0-|-0 colliding wind binary. We showed 
that the X-ray flux emitted by HD 152248 presents a clear 
enhancement of about 50 per cent culminating slightly 
after apastron passage and that the relative amplitude of 
this increase is also dependent on the considered energy 
range. Comparison with previous ROSAT PSPC data 
provides a good support for a phase-locked behaviour of the 
observed modulations that further remains stable on time 
scales of at least 10 years. The XMM EPIC spectra are 
relatively soft and present their maximum peak between 
0.8 and 0.9 keV. They are reasonably well described by a 
two-temperature mekal model in the range [0.5 -2.5 keV]. 
The temperature of the low-energy component is located 
between 0.26 and 0.31 keV while the second component 
has a temperature ranging from 0.62 to 0.71 keV (see 
Table 5). A third component of higher energy is clearly 
seen in several spectra, though its nature could not be 
unambiguously established. It could indeed correspond to 
a power-law tail characterized by a photon index F ~ 3.7 
or, though with a slightly lower confidence, to a thermal 
component with a temperature of a few keV. If the existence 
of the non-thermal component is indeed confirmed, then its 
spectral index is much larger than the 1.5 value expected for 
a strong shock. Indeed, Chen & White (1991) expect such a 
component to be produced by inverse Compton scattering 
radiation emitted by a population of relativistic electrons 
that would have been accelerated in shocks through the 
first-order Fermi process. This result is similar to what 
has been observed for other early type O stars such as 
9 Sgr and HD 93403 (Rauw et al. 2002a,b) and the problem 
definitely deserves further attention. Combining the data 
from all six pointings however reveals that the Fe line at 
6.7 keV is most probably present. This implies that at least 
part of the high energy tail is likely to have a thermal origin. 

We also performed numerical hydrodynamical simu- 
lations that closely reproduce the geometrical and orbital 
configuration of the HD 152248 system at the six phases 
corresponding to the different XMM-Newton pointings. The 
resulting hydrodynamical maps reveal that the interaction 
region is highly unstable and, as we progress away from 
the binary axis, it develops a filamentary structure that 
alternates thin but dense shells with larger low density 
regions. The computed averaged intrinsic flux produced 
by such a configuration is in excellent agreement with the 
analytical estimate for the X-ray emission resulting from 
the interaction region of two identical winds. For a given 
configuration, the absorbed X-ray luminosity predicted by 
the model further present clear variations on time scales of 
a couple of hours. Finally, the predicted luminosities show a 
clear enhancement by about one order of magnitude around 
apastron. Adding the expected intrinsic contribution of the 
two O-type stars of the system, the total X-ray luminosity 
reasonably well reproduces both the observed luminosities 
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and the amplitude and phase of their variation. Due 
to the 2-D nature of the simulations and the fact that 
we neglect the orbital motion, the predicted light curve 
is symmetric around apastron and could therefore not 
properly render the observed decrease of the flux beyond 
periastron phase. Part of the answer to this problem 
may probably be found in the deflection of the interaction 
region due to the orbital motion and subsequent absorption. 

The energy distribution predicted by the hydrody- 
namical simulations are peaked around 3keV and are 
therefore much harder than the observed EPIC spectra. 
The models undcr-predict the flux at low-energy, which 
is most probably due to the fact that the X-ray contri- 
bution of the two O-stars has been neglected. Similarly, 
the hard-energy tail seems to be overestimated in the 
model prediction, though the observed energy distribution 
still remains within the boundaries of the model. Wo are 
conscious that our simulations rest on several oversimplified 
hypotheses such as constant wind velocities and frozen or- 
bital configurations, whose efiiects may influence the output 
of our models. We are however confident that our results 
provide a good approximation to the 'reality'. Another 
limitation that we are facing is the extremely thin inter- 
action region that could not be resolved with the adopted 
numerical mesh. One significant improvement for hydro- 
dynamical simulations of highly radiative colliding wind 
interactions could therefore be the use of adaptative meshes. 

Using difi^erent techniques, we also investigated the 
X-ray variability of the source on short-time scales. Linear 
trends are definitely present in several data sets and helped 
to constrain the first order derivative of the flux. We also 
applied the Detrended Fluctuation Analysis, a method 
that has proved its merits in other domains of sciences. 
Its application to the EPIC data suggests that long-range 
correlations might be present in the series formed by 
the photons arrival times, especially while restricting 
the analysis to the hard energy tail of the data. This 
however does not seem to be directly related to possible 
turbulent processes occurring within the interaction region. 
Short-time variability with time scales of a couple of hours 
as seen in the predictions of the simulations could finally 
not be brought into light. On one hand, it is possible that 
these variations are diluted by the additional flux from the 
two stars or perturbed by some physical phenomenon not 
accounted for in the simulations, such as the orbital motion. 
On the other hand, the predicted variability might be an 
artefact resulting from the 2-D nature of the simulations in 
which the symmetry might lead to an overestimate of the 
variability. Again 3-D simulations are required to clarify 
this question. 

Finally, wo emphasize that this campaign on the 
HD 152248 massive binary is one of the most complete 
phase-resolved X-ray studies of a colliding wind O-l-O 
binary ever performed. It further outlines the insight in 
the physics of the phenomenon that can be gained from 
a detailed comparison of the observations with hydrody- 
namical models. It also emphasizes that more detailed 
hydrodynamical models, possibly using 3D configuration 
and accounting for the orbital motion, are needed to 



achieve a better understanding of the local physics that 
governs the X-ray emission of colliding wind binaries. The 
need to extend this kind of study to a growing number 
of CWB systems, either 0-1-0 or WR-fO binaries, is 
obvious to achieve a detailed understanding both of the 
hot star winds and of their interaction within CWB systems. 
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APPENDIX A: THE PROBABILITY OF 
VARIABILITY TEST REVISITED 

Al Preibisch & Zinnecker pov test 

The probability of variability test suggested by Preibisch & 
Zinnecker (2002) basically gives the confidence level - called 



Table Al. Example of the dependence of the pov estimator with 
the number of observed intervals nint . First column gives the val- 
ues of riint ■ Second to fourth columns provide the pov confidence 
level for a maximum observed counts A^max respectively of 5, 6 

and 8 as computed following our revised version with < TV >= 2 
(see Sect. A2). Preibisch & Zinnecker would have obtained values 
corresponding to nint = 1 for each case. 



nint probability of variability 





-^max — 5 


A''max = 6 


A^^max — 8 


1 


0.9473 


0.9834 


0.9989 


2 


0.8975 


0.9671 


0.9978 


5 


0.7630 


0.9197 


0.9945 


10 


0.5822 


0.8458 


0.9891 


20 


0.3390 


0.7154 


0.9783 


100 


0.0045 


0.1882 


0.8961 


1000 


0.0000 


0.0000 


0.3338 



the probability of variability (pov) - to reject the null hypoth- 
esis of a constant count rate for the considered source. To do 
so, given the mean count rate, they first estimate the mean 
time interval At during which two photons are expected to 
be detected by the instrument. Then, using a sliding win- 
dow, they determine the maximum number of counts iVmax 
received during such a period At. Finally, they estimate the 
probability p that a poissonian distribution of mean 2 gives 
at least A''max counts in a At interval. This is simply given 
by the formula ; 

p=l-g = l- £ e-'^ (Al) 

The probability of variability given by pov = 1— p estimates 
whether such a high number of counts Nuiax as observed is 
rather improbable, which is expressed by a large value of 
the pov and suggests a significant variability in the data. 

However, in their method, Preibisch & Zinnecker do 
not account for the fact that, actually, they are not per- 
forming a single experiment (i.e. observing a single interval) 
but rather riint of them, where the value for nint can be 
quite large (a few thousands in the case of our present 
observations). Indeed, repeating a counting experiment 
increases the probability to observe a number of counts that 
seems to significantly deviate from the expected values but 
still is in fair agreement with the statistics theory. Preibisch 
& Zinnecker also failed to consider a more subtle aspect: 
the fact that their "sliding window" approach implies 
that the observed intervals are overlapping and thus that 
the counting of photons in the different windows are not 
independent. In its current form, the pov confidence level 
used by Preibisch & Zinnecker is biased towards larger 
values (see Table Al for an example) and clearly favours 
the rejection of the null hypothesis of non variability, which 
might lead to a dramatic increase of type I error occurrences. 



A2 The pov revisited 

Let At be the time period during which the expected number 
of counts is<iV> (<Ar>=2in Preibisch & Zinnecker 
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original method). Let nint be the number of non overlapping 
intervals of size At among which we distributed our data. 
Finally, let A'max be the maximum number of counts found in 
a single interval among the riint ones. Therefore, as the nint 
countings might reasonably be considered as independent 
and neglecting the fact that At has been estimated from 
the time series, the probability to find by chance (i.e. due to 
statistical fluctuations) at least iVmax counts in one or more 
intervals among the rii^t observed ones is now given by: 

while the probability of variability, i.e. the probability that 

such a high rmmbcr of counts A'max is not duo to statistical 
fluctuations, is given by pov = 1 — P. As of course g < 1, 
P increases with the number of intervals nint and the pov 
might thus be drastically reduced (see Table Al). 

In the regard of the present considerations, the pov 
levels found by Preibisch & Zinnecker (2002) in their 
Section 7 on the variability of the studied objects are 
therefore overestimated or, at least, could not be compared 
to a unique rejection criterion as larger values of their pov 
confidence level are naturally expected for data displaying a 
larger number of counts. This is especially true concerning 
the objects with the higher number of counts (N-045-02 and 
L-312) for which the large values of their pov confidence 
level found {pov = 0.95) are most probably due to statistical 
fluctuations rather than reflecting any intrinsic variability 
of the sources. 



